/**************************************************************************************
            UNIVERSIDADE FEDERAL DE PERNAMBUCO

Author:     Rogerio Soares da Silva
Date:       May-2006
Licence:    This program is free software; you can redistribute it and/or modify it.   
Warranty:   None! Use it at your own risk!
Visit:      www.padmec.org

 **************************************************************************************/
#include "manager.h"   


class IdPair
{
public:
	IdPair(int id0,int id1):_id0(id0),_id1(id1){}
	~IdPair(){}
	int getId0() const { return _id0; }
	int getId1() const { return _id1; }

private:
	int _id0;
	int _id1;
};


// local class for unify Cij's
class EdgeInfo
{
public:

	enum xyzLists {xlist, ylist, zlist};
	typedef list<IdPair>::iterator IDIter;
	typedef list<double>::iterator CijIter;
	typedef map<int, int >::const_iterator MCI; // MCI - map const iterator

	EdgeInfo() {}

	~EdgeInfo() {}

	IDIter IdPair_beginIterator() { return listEdgeIds.begin(); }
	IDIter IdPair_endIterator() { return listEdgeIds.end(); }

	void printEdgeIds()
	{
		CijIter Cijx = xList.begin();
		CijIter Cijy = yList.begin();
		CijIter Cijz = zList.begin();
		for (IDIter iter = listEdgeIds.begin(); iter != listEdgeIds.end(); iter++,Cijx++,Cijy++,Cijz++)
			printf("[%d]aresta: %d %d :: %f %f %f\n",P_pid(),iter->getId0()+1, iter->getId1()+1,*Cijx,*Cijy,*Cijz);
	}

	CijIter Cij_beginIterator(xyzLists whichone)
	{
		whichList = whichone;
		return Cij_beginIterator();
	}

	CijIter Cij_beginIterator()
	{
		switch (whichList)
		{
		case xlist:  return xList.begin(); break;
		case ylist:  return yList.begin(); break;
		case zlist:  return zList.begin(); break;
		}
	}

	CijIter Cij_endIterator(xyzLists whichone)
	{
		whichList = whichone;
		return Cij_endIterator();
	}

	CijIter Cij_endIterator()
	{
		switch (whichList)
		{
		case xlist:  xList.end(); break;
		case ylist:  yList.end(); break;
		case zlist:  zList.end(); break;
		}
	}

	void printListVals()
	{
		CijIter Cijx = xList.begin();
		CijIter Cijy = yList.begin();
		CijIter Cijz = zList.begin();
		for (; Cijx != xList.end(); Cijx++,Cijy++,Cijz++)  printf("[%d] Cij: %f %f %f\n",P_pid(),*Cijx,*Cijy,*Cijz);
	}
	int numEdges() const {  return listEdgeIds.size();  }

	void setWhichList(xyzLists wl) { whichList = wl; }

	void setEdgeIds(IdPair ids)
	{
		listEdgeIds.push_back(ids);
	}
	const list<IdPair>& getRemoteEdgeIds() { return  listEdgeIds; }

	const list<double>& getRemoteEdgeCijList()  { return getRemoteEdgeCijList(whichList); }

	void setRemoteEdgeCij(const double &val)
	{
		setRemoteEdgeCij(whichList,val);
	}
	void setRemoteEdgeCij(xyzLists whichone, const double &val)
	{
		switch (whichone)
		{
		case xlist:  xList.push_back(val); break;
		case ylist:  yList.push_back(val); break;
		case zlist:  zList.push_back(val); break;
		}
	}

	int sizeRemoteEdgeCij(xyzLists whichone)
	{
		switch (whichone)
		{
		case xlist:  return xList.size(); break;
		case ylist:  return yList.size(); break;
		case zlist:  return zList.size(); break;
		}
	}

	const list<double>& getRemoteEdgeCijList(xyzLists whichone)
	{
		switch (whichone)
		{
		case xlist:  return xList; break;
		case ylist:  return yList; break;
		case zlist:  return zList; break;
		}
	}

	void updateRemoteEdgeCijList(CijIter iter, double val)
	{
		switch (whichList)
		{
		case xlist:  xList.insert(iter, val); break;
		case ylist:  yList.insert(iter, val); break;
		case zlist:  zList.insert(iter, val); break;
		}
	}

	// after everything is done, clean the room
	void deleteEdgeIds() { listEdgeIds.clear(); }

	// after update each coordenate of Cij delete it
	void deleteCij_list()
	{
		deleteCij_list(whichList);
	}
	void deleteCij_list(xyzLists whichone)
	{
		switch (whichone)
		{
		case xlist:  xList.clear(); break;
		case ylist:  yList.clear(); break;
		case zlist:  zList.clear(); break;
		}
	}

	//mapRows() must be called BEFORE use numRowstoImport() and rowstoImport()
	void mapRows()
	{
		set<int> sTemp;		
		for ( IDIter iter = listEdgeIds.begin(); iter != listEdgeIds.end(); iter++)
			sTemp.insert( iter->getId0() );

		int i = 0;
		for( set<int>::iterator sIter = sTemp.begin(); sIter != sTemp.end(); i++, sIter++)
			mRows.insert(make_pair(*sIter, i) );
		sTemp.clear();
	}

	void showMapping()
	{
		MCI miter = mRows.begin();
		for (; miter !=mRows.end(); miter++)
			printf("[%d] ido=%d -> %d\n",P_pid(),miter->first, miter->second);
	}

	int numRowstoImport() const { return mRows.size(); }

	void rowstoImport(int *rows) const
	{
		MCI mapIter = mRows.begin();
		for ( int i=0; mapIter != mRows.end(); mapIter++, i++)
			rows[i] = mapIter->first;
	}

	// which row must be read from parallel matrix
	int findMappedRow(const int &row)
	{
		MCI mapIter = mRows.find(row);
		return mapIter->second;
	}

	void deleteMapRows() { mRows.clear(); }

private:

	xyzLists whichList;               // enum to specify which list of coordenate must be returned
	list<IdPair> listEdgeIds;     // list for edges' ids
	// list to store Cij coordenates
	list<double> xList;
	list<double> yList;
	list<double> zList;

	map<int,int> mRows;   // mapping to chose which rows must be import from parallel matrix

};

#ifdef _HAVE_PETSC_
int unifyCijAmongProcessors(EdgeInfo* edgeinfo, pMesh theMesh);
// ===================================================================================


// In parallel Cij is calculated partialy on each partition, but they need the same value
// as they were calculated in serial. Cij(rank p, dom k) = Cij(rank q, dom k)
//
// Cij(rank p, dom k) = Cij(rank p, dom k) + Cij(rank q, dom k)
// Cij(rank q, dom k) = Cij(rank q, dom k) + Cij(rank p, dom k)
//
// Only edges on partion boundary are treated here. Do nothing in serial.
// =========================================================================================
// IN:
// void
// OUT:
// void
int unifyCijAmongProcessors(pMesh theMesh, const set<int> &setOfDomains)
{
	// GO AHEAD ONLY IF THERE ARE MORE THAN ONE PROCESSOR

	if (P_size()==1) return 1;
	printf("PARALLEL: %s\n",__FILE__);

	EdgeInfo* edgeinfo;
	char *Cij_xyz[3] = {"Cx-","Cy-","Cz-"};

	int i;
	pEntity edge;
	vector<pVertex> vertex;

	set<int>::iterator SIter;
	for (SIter = setOfDomains.begin(); SIter != setOfDomains.end(); SIter++)
	{
		const int dom = *SIter;	// domain flag		
		edgeinfo = new EdgeInfo;	// used to store edge node's IDs and Cij coordenates

		// loop over all edges: get only those with remote copies
		// get edge node's IDs and Cij
		EIter eit = M_edgeIter(theMesh);
		while (edge  = EIter_next(eit))
		{
			if ( M_numRemoteCopies(theMesh,edge) > 0 )
			{
				M_GetVertices(edge,vertex);
				int id0 = EN_id(vertex[0])-1;
				int id1 = EN_id(vertex[1])-1;
				if (id1 < id0) swap(id0,id1);
				vertex.clear();

				// GET IDS
				IdPair idpair(id0,id1);
				edgeinfo->setEdgeIds(idpair);

				// GET Cij
				double val[3] = {.0,.0,.0};
				char CijString[256];				
				for (i=0; i<3; i++)
				{
					sprintf(CijString,"%s%d",Cij_xyz[i],dom);						
					EN_getDataDbl(edge,MD_lookupMeshDataId(CijString),&val[i]);
				}

				edgeinfo->setRemoteEdgeCij(edgeinfo->xlist,val[0]);
				edgeinfo->setRemoteEdgeCij(edgeinfo->ylist,val[1]);
				edgeinfo->setRemoteEdgeCij(edgeinfo->zlist,val[2]);
			}
		}
		edgeinfo->mapRows();    // não comente esta linha!!!!

		//edgeinfo->showMapping();
		//edgeinfo->printEdgeIds();

		cout<<__LINE__<<endl;
		// UNIFY CIJ FOR EDGES WITH REMOTE COPIES
		// ---------------------------------------------------------------------
		// unify Cij_x

		edgeinfo->setWhichList(edgeinfo->xlist);
		unifyCijAmongProcessors(edgeinfo,theMesh);

		// unify Cij_y
		edgeinfo->setWhichList(edgeinfo->ylist);
		unifyCijAmongProcessors(edgeinfo,theMesh);

		// unify Cij_z
		edgeinfo->setWhichList(edgeinfo->zlist);
		unifyCijAmongProcessors(edgeinfo,theMesh);


		// update Cij on mesh data structure
		// ---------------------------------------------------------------------
		list<double>::iterator Cij_xIter = edgeinfo->Cij_beginIterator(edgeinfo->xlist);
		list<double>::iterator Cij_yIter = edgeinfo->Cij_beginIterator(edgeinfo->ylist);
		list<double>::iterator Cij_zIter = edgeinfo->Cij_beginIterator(edgeinfo->zlist);

		// loop over all edges once more to update only those with remote copies
		eit = M_edgeIter(theMesh);
		while (edge  = EIter_next(eit))
		{
			if (M_numRemoteCopies(theMesh,edge) > 0) // only edge with remote copies
			{
				char CijString[256];

				sprintf(CijString,"%s%d",Cij_xyz[0],dom);
				EN_attachDataDbl(edge,MD_lookupMeshDataId(CijString),*Cij_xIter); Cij_xIter++;

				sprintf(CijString,"%s%d",Cij_xyz[1],dom);
				EN_attachDataDbl(edge,MD_lookupMeshDataId(CijString),*Cij_yIter); Cij_yIter++;

				sprintf(CijString,"%s%d",Cij_xyz[2],dom);
				EN_attachDataDbl(edge,MD_lookupMeshDataId(CijString),*Cij_zIter); Cij_zIter++;
			}
		}
		EIter_delete(eit);

		//cout<<__LINE__<<endl;

		// free memory
		edgeinfo->deleteEdgeIds();
		edgeinfo->deleteMapRows();
		edgeinfo->deleteCij_list(edgeinfo->xlist);
		edgeinfo->deleteCij_list(edgeinfo->ylist);
		edgeinfo->deleteCij_list(edgeinfo->zlist);
		delete edgeinfo; edgeinfo = 0;

	} // END OF LOOP OVER DOMAINS

	MPI_Barrier(MPI_COMM_WORLD);
}

// Calculates geometric coefficient Dij and assign it to faces (IT WILL BE MODIFIED SOON)
// =========================================================================================
// IN:
// edgeinfo - EdgeInfo object with a list of remote edges that must be unified
// OUT:
// void
int unifyCijAmongProcessors(EdgeInfo* edgeinfo, pMesh theMesh)
{   
	PetscErrorCode ierr;
	double val;
	int i;
	int numNodesGlobal = M_getMaxVId(theMesh);	

	int d_nz = (int)(numNodesGlobal)/P_size();
	int o_nz = (int)(numNodesGlobal)*(P_size()-1)/P_size();

	// matrix to store Cij coord for all processes
	Mat exchangeCij; // matrix to store nodes from all processes
	ierr = MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE,PETSC_DECIDE,numNodesGlobal-1,numNodesGlobal,0,PETSC_NULL, 0,PETSC_NULL, &exchangeCij); CHKERRQ(ierr);

	// transfer all Cij coordenates to a parallel matrix. duplicated Cij are added      
	// maybe CijCoords can be inserted into exchangeCij without a loop.
	const list<IdPair> idpairList = edgeinfo->getRemoteEdgeIds();

	list<IdPair>::iterator Id_Iter = edgeinfo->IdPair_beginIterator();
	list<double>::iterator Cij_Iter = edgeinfo->Cij_beginIterator();

	// populate exchangeCij matrix
	for (; Id_Iter != edgeinfo->IdPair_endIterator(); Id_Iter++, Cij_Iter++)
	{
		ierr = MatSetValue(exchangeCij,Id_Iter->getId0(),Id_Iter->getId1(),*Cij_Iter,ADD_VALUES); CHKERRQ(ierr);
	}

	ierr = MatAssemblyBegin(exchangeCij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
	ierr = MatAssemblyEnd(exchangeCij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
	//ierr = MatView(exchangeCij,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

	// now, all nodes, regardless their location(processor), are updated.
	Mat unifiedCij;     // matrix with all requested nodes

	PetscInt nrows = edgeinfo->numRowstoImport();           // how rows rank p must require from parallel matrix
	//  printf("[%d] number of rows: %d\n",P_pid(),nrows);                                            // each process can only read owned rows

	PetscInt ncols = numNodesGlobal;

	int* rows = new int[nrows];                                     // rows to be imported
	edgeinfo->rowstoImport(rows);

	//    for(i=0; i<nrows; i++) printf("[%d] rows: %d\n",P_pid(),rows[i]);

	int* cols = new int[ncols];
	for (i=0; i<numNodesGlobal; i++) cols[i] = i;

	// coefficients will now be accessable locally
	MatGetSubMatrixRaw(exchangeCij, nrows, rows, ncols, cols, PETSC_DECIDE,	MAT_INITIAL_MATRIX, &unifiedCij);  CHKERRQ(ierr);

	ierr = MatDestroy(exchangeCij); CHKERRQ(ierr);
	delete[] rows;
	delete[] cols;
	//    ierr = MatView(unifiedCij,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

	PetscInt m, n, row, col;
	ierr = MatGetOwnershipRange(unifiedCij,&m,&n);  // range of rows from parallel matrix owned to rank p

	// loop to update Cij coefficients list
	//printf("rank %d\tM = %d N = %d\n",P_pid(),m,n);
	edgeinfo->deleteCij_list();
	for (Id_Iter = edgeinfo->IdPair_beginIterator(); Id_Iter != edgeinfo->IdPair_endIterator(); Id_Iter++, Cij_Iter++)
	{   
		row = edgeinfo->findMappedRow( Id_Iter->getId0() ) + m;
		col = Id_Iter->getId1();
		//       printf("[%d] aresta: [%d(mapped: %d), %d]\n",P_pid(),Id_Iter->getId0(),row,Id_Iter->getId1());
		ierr = MatGetValues(unifiedCij,1,&row,1,&col,&val);CHKERRQ(ierr);
		//         printf("[%d] edge[%d - %d] = %f\n",P_pid(),Id_Iter->getId0(),col,val);
		//edgeinfo->updateRemoteEdgeCijList(Cij_Iter,val);
		edgeinfo->setRemoteEdgeCij(val);
	}
	ierr = MatDestroy(unifiedCij); CHKERRQ(ierr);
}
#endif //-D_HAVE_PETSC_
